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Abstract. As part of our development of a computer code to perform 3D 
'constrained evolution' of Einstein's equations in 3+1 form, we discuss issues 
regarding the efficient solution of elliptic equations on domains containing holes 
(i.e., excised regions), via the multigrid method. We consider as a test case the 
Poisson equation with a nonlinear term added, as a means of illustrating the 
principles involved, and move to a "real world" 3-dimensional problem which 
is the solution of the conformally flat Hamiltonian constraint with Dirichlet 
and Robin boundary conditions. Using our vertex-centered multigrid code, we 
demonstrate globally second-order-accurate solutions of elliptic equations over 
domains containing holes, in two and three spatial dimensions. Keys to the 
success of this method are the choice of the restriction operator near the holes 
and definition of the location of the inner boundary. In some cases (e.g. two 
holes in two dimensions), more and more smoothing may be required as the mesh 
spacing decreases to zero; however for the resolutions currently of interest to 
many numerical relativists, it is feasible to maintain second order convergence by 
concentrating smoothing (spatially) where it is needed most. This paper, and our 
publicly available source code, are intended to serve as semi-pedagogical guides 
for those who may wish to implement similar schemes. 



1. Introduction 

Solving Einstein's equation in 3+1 form jTl |2] requires that a set of elliptic (or 
quasi-elliptic) constraint equations be satisfied at all times. The Bianchi identities 
ensure that, given a set of initial data which analytically satisfy the constraints, the 
subsequent analytically evolved variables will also satisfy the constraints. In numerical 
solutions of Einstein's equations, however, the constraints are not preserved exactly. 
Thus errors will arise as the simulation proceeds, and the extent to which the numerical 
solutions actually reflect the true solutions of the analytic equations is still an open 
area of research (e.g., E El EHZ]). Apart from the issue of the accuracy of the 
solutions obtained, there is also the problem of numerical (in) stability, which seems 
to be related to the lack of preservation of the constraints. 

It has been observed for lower-dimensional computational relativity simulations 
|H1 El d| that the use of "constrained evolution" schemes, in which the constraint 
equations are satisfied (to some fixed accuracy) at all timesteps, can offer better 
stability behavior. Such schemes have not yet received adequate attention in 3D 
simulations simply because of the costly nature of solving the constraints at every 
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timestep. We are among a set of researchers (e.g., ^2 El El) interested in exploring 
the benefits of constrained evolution for 3D applications. 

Our desire to solve the constraints at each timestep motivates us to develop a 
fast elliptic solver with which to enforce the constraints, and we turn our attention to 
one of the most popular methods currently in use: the multigrid method The 
multigrid method is particularly attractive because it is an optimal method, i.e. it 
requires only 0{N) operations, where A'^ is the number of unknowns. Furthermore, 
the multigrid method is not especially difficult to implement. (Although multidomain 
pseudospectral collocation methods ^1 El also offer very efficient solutions of elliptic 
problems, we focus on multigrid because of the relative ease with which one can develop 
a multigrid code. It is not our intent in this paper to present a comparison between 
multigrid methods and other numerical methods.) In the course of developing our own 
implementation of the multigrid method containing holes (i.e., excised regions), it came 
to our attention that some researchers consider it impossible, or at least infeasible in 
practice, to obtain generic results which have everywhere the same order of accuracy 
as the finite difference scheme employed, for domains containing holes, particularly in 
three dimensions. (Specifically, that, using second-order-accurate difference stencils, 
one could not generically obtain second-order-accurate multigrid solutions if part of 
the domain was excised.) 

Since we expect to use exact, analytic data as an inner Dirichlet boundary 
condition (at the surface of the hole), we believe it is not necessary to employ 
'sophisticated' schemes such as deferred correction ^21, which can be used to achieve 
higher order accuracy for various boundary conditions. We are already aware that 
one can fairly easily obtain second order convergence using square (in 2D) or cubical 
(in 3D) excision regions in which the physical domain of the hole is the same on all 
multigrid levels. However, we wish to apply our multigrid code in situations where the 
shapes of the holes are more general (e.g., spherical holes). Thus we are motivated to 
find a simple, robust method to handle the physical situations of interest to us. 

If the multigrid solver is to be used to provide initial data for a numerical evolution 
code, the error of the initial solution need only be below the truncation error of the 
finite difference scheme used for evolution ^2]. Thus our pursuit of specifically second 
order convergence for the multigrid scheme might be regarded as unnecessary. Our 
motivation is partly to develop a simple algorithm which nevertheless offers fast (better 
than first order) convergence, and partly because, given our exact Dirichlet conditions 
for the inner boundary, we simply expect that it should not be difficult to obtain 
second order convergence from the multigrid solver. 

Although it is likely the case that some of the techniques presented here are 
known to some specialists, we believe it is worthwhile to popularize the straightforward 
methods we have used. 

If we were applying a numerical method to a new physical system, we would want 
to provide a full discussion of mathematical foundations, for issues such as existence 
and uniqueness of the solutions. Since the use of finite difference techniques to solve the 
constraint equations is by now a standard activity in relativity, and the properties of 
the constraint equations are well known (for the Euclidean, constant-mean-curvature 
background and for the signs of the solutions we use in this paper f^l 1201 1^^- 
we will focus our attention on the benefits of our multigrid implementation: second 
order accuracy near the hole via a simple method, and the efficiency improvements 
afforded by concentrating smoothing near the hole and performing more smoothing 
sweeps during the pre-Coarse-Grid-Correction (CGC) phase than during the post- 
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CGC phase. The vaUdity of our results wiU be estabhshed by demonstrating second 
order convergence to known exact solutions as well as second order behaviour in the 
truncation error. 

The organization of this paper is as follows: We begin with an overview of 
the multigrid method. We then discuss some results obtained in two dimensions, 
which we then extend toward a more difhcult equation in three dimensions, eventually 
demonstrating the method by solving the Hamiltonian constraint equation for a black 
hole (with additional matter) on a Euclidean background. Finally we share some ideas 
for more general applications, and we end with a summary of the main five 'tips' we 
have for others wishing to implement multigrid methods on domains with holes. 

2. Overview of Multigrid 

The multigrid method, first introduced by Brandt has received considerable 

attention in the literature, and is the subject of numerous articles, conferences, reviews 
and books (e.g., |22II23[I^I25[I^I27| ). However, its application on domains with holes 
— or on domains with "irregular boundaries" in general — has received only modest 
attention. The excellent works of Johansen and Colella "SS* and Udaykumar et al. 23 
for cell-centered multigrid are notable exceptions. The "BAM" code [HOI is a notable 
introduction of the multigrid method in numerical relativity, and even features the 
ability to handle domains with multiples holes. (It does not employ the advances we 
describe in this paper, notably our method of obtaining second order accuracy near 
the hole.) Despite this significant result, multigrid methods in numerical relativity 
remain as tools of only a few specialists. We thus thought it worthwhile to share our 
experiences and methods with vertex-centered multigrid with the numerical relativity 
community. 

A simple algorithm to solve discretized elliptic equations is Gauss-Seidel 
relaxation [31] ■ This method works very poorly at high resolutions because it fails 
to operate efficiently on long-wavelength components of the error. However, it is 
extremely effective at eliminating short-wavelength components of the error, or in 
other words, at "smoothing" the error (i.e., the residual, see below). The multigrid 
scheme is essentially a clever means of eliminating successive wavelength-components 
of the error via the use of relaxation at multiple spatial scales. 

Here we give a very brief overview of the multigrid method, following the notes by 
Choptuik 32 . (Introductions to multigrid applications in numerical relativity are also 
found in Choptuik and Unruh |33| and Brandt |34j.l We want to solve a continuum 
differential equation Cu — /, where £ is a differential operator, / is some right hand 
side, and u is the solution we wish to obtain. We discretize this differential equation 
into a difference equation on some grid (or lattice) with uniform spacing h: 

= f\ (1) 

where is the exact solution of this discrete equation, and lim/j^o w'' = u. Rather 
than attempting to solve ^ directly via the costly operation of matrix inversion, we 
apply an iterative solution method. At any step in our iteration, we will have only an 
approximate solution 2± u'', such that 

^h^h _ jh ^ ^2) 

where is some small quantity called the residual. The variables and ■u'* are 
related by 

u''^u''+v'', (3) 
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where is some correction term. In this iterative algorithm, we start with some guess 
^oid ^'^"^ ^'^ bring it closer to by applying some (approximate) correction: 

"ncw:="old+^' (4) 

The remainder of the scheme involves making a choice of what to use for v^. A simple 
choice, called the Linear Correction Scheme (LCS), can be used for hnear operators 
C^. This method can then be extended to work with nonlinear operators via the Full 
Approximation Storage (FAS) method, below. 

2.1. Linear Correction Scheme (LCS) 

For linear £, we begin by substituting Q into (QJ: 

C''\u^ + v^) ^ . (5) 
Then we use the linearity of C: 

C^u^ + C^v^ = f^, (6) 

= -r\ (7) 
Consider a coarser- grid version of 10, in which we use a mesh spacing of 2/i: 

Now we come to one of the 'tricks' of the multigrid method: Instead of the coarser 
residual in equation JHJ, use 

f2h^j2h~h^ (9) 

where is a restriction operator, which maps values from the fine grid to the coarse 
grid via some weighted averaging operation. So then we have 

C^h~2h ^_j2h~h^ (10) 

(Note the tilde on u^'*, denoting that this is not the exact v^^ , because we use a 
residual restricted from the fine grid.) Equation H1U|) can be solved "exactly" for ifl^ 
because this is inexpensive to do on the coarse grid. 

A second 'trick' is that, for our correction term in our numerical update scheme 
(@J, we do not directly solve equation (|7I), but instead we use u'' approximated from 
the coarser (2ft,) correction iP'^: 

v''=lLv'\ (11) 

where /j/j is an interpolation or prolongation operator, which maps values from the 
coarse grid to the fine grid via some interpolation operation. Thus we use C^'* as a 
coarse grid correction (CGC) to w'': 

"now •= "old + ^211 ^"^^^ (12) 

In performing the restriction f^'' = if^ , we are assuming that is sufficiently 
smooth to be sensibly represented on the coarse grid (e.g., without aliasing effects). 
This implies that vP' also needs to be smooth for this restriction to produce meaningful 
results. Therefore, before each restriction operation, we apply a series of "smoothing 
sweeps" to vP" in an effort to smooth the residual f'', using the efficient smoothing 
algorithm of Gauss-Seidel relaxation. 
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2.2. Full Approximation Storage (FAS) method 

For nonlinear operators £'', we must modify the algorithm outlined above. Our 
implementation relies on the so-called "alternative" description of the FAS algorithm, 
which involves a notion of the truncation error t^^ , defined on the coarse grid by 

t'"' =C^^u- f^, (13) 

where u is the exact solution to the continuum equation. The function t'^^ can be 
regarded as a correction term which makes the finite difference equation produce the 
continuum solution: 

C^hu^ + r^h^ (14) 

For general problems, the continuum solution u and truncation error r^'' may not be 
available, however we may use an approximation to t'^^ which is the relative truncation 
error between the coarse and fine grids, r^'*, given by 

tI^-C^'^II^u'^-II'^C'^u^. (15) 

Using this in H14|l gives 

r2h„,2h ^ r2h . 2h 

L u ~ / +Tf^ , (Id) 

and it is this equation which we solve on the coarse grid. 

To obtain the CGC to a fine grid solution, we substitute H15|l into H16|) . obtaining 

^2h^2h _ ^2hj2h ~h ^ j2h ^jh _ jrh~h^^ (^7) 

where the term l'^'^ on the right is introduced for f"^^, in analogy with equation 
(|SJ). By analogy to H1U|) in the LCS scheme, the term we should use for the CGC is 
(the part on the left hand side of the previous equation, without the £'s): 

~2h ^^2h_j2h^h^ (18) 

Thus we arrive an update scheme of 

uL^:^u'M + l2Hiu"'-ll''u'^)- (19) 

Note that u^^ is the exact solution to 1)16(1 . and can be obtained with little effort due 
to the low resolution on the coarse grid. (If there are more than two multigrid levels 
involved in the solution, then it^'* should be the best approximation u^'^ obtained on 
the coarser grid.) 



2.3. V- Cycles and the Full Multigrid Algorithm 

Instead of only using two grids as described above, one could find the coarse grid 
correction (CGC) to a fine-grid-problem by solving for a CGC from an even coarser 
grid, i.e., obtain -u^'' ~ u^'* by finding u'^^ (in which case the corresponding right hand 
side is /'*''+T2^', where = ^2^'[/^''+t^''] ). One can imagine a hierarchy of multiple 
such grids, in which coarser grids provide CGCs for finer grids. The solution algorithm 
will then take the form of a V- cycle., in which we start with an initial guess on the fine 
grid, at multigrid level /max- Then we perform some number of smoothing sweeps and 
restrict the data to a coarser grid. We continue smoothing and restricting to coarser 
grids until we arrive at a grid coarse enough to solve the coarse grid equation ((TE|l 
'exactly' (i.e., to machine precision), at minimal computational cost. This coarsest 
grid is at level Inun- We then prolongate this solution to finer grids by performing a 
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series of coarse-grid corrections, with perhaps additional smoothing operations being 
performed before moving to each finer grid. 

In addition, before starting a V-cycle from the finest grid, we can use an initial 
guess obtained from a prior solution at a coarser resolution. Doing this for each grid 
level results in the Full Multigrid Algorithm (FMA). We outline the FMA as follows, 
(using a notation where superscripts refer to multigrid levels, with /min the coarsest 
level, and ^max the finest): 



u = Mo 
Solve L'—u' 

do I = ^min + 



1 to 



do TO = ^ to L 



'exactly' 



1 



Smooth (solve via Gauss-Seidel) at level to 



end do 



I" 



-1 



II Initial guess, e.g. mq = 
// or uo = 1 



// Initial guess for fine grid 
// Begin V-cycle 

// Pre-CGC smooths 

// Restrict to coarser grid 

// Restrict to coarser grid 



Solve L'""" u' 



/'■"■" 'exactly' 



// Solve on coarsest grid 



do TO = ^ 
,1 



mm 



1 to L 



:= + C_i [u- 
Smooth at level to 
end do 



// Apply CGC 

// Post-CGC smooths 

// End V-cycle 
// End FMA 



end do 

Thus we see that on all grids except the coarsest grid, we only smooth the error, 
and we solve the difference equation exactly only on the coarsest grid. 



3. Solution of a Nonlinear Poisson Equation 

Since we are interested in ultimately solving for the constraint equations in general 
relativity, we consider the Haniiltonian constraint in a conformally flat background 
geometry. This yields the the Poisson equation with a pair of nonlinear terms added: 

^2 

V^u(a;, y, z) - K^u^(x, y, z) + = j(x, y, z), (20) 

u'{x,y,z) 

where is the usual Laplacian in Euclidean space. and are arbitrary positive 
real constants related to the rate of expansion for the 3-space for which is 120() is the 
constraint equation. f{x, y, z) is related to the energy density in the 3-space, and can 
be chosen such that the resulting u{x, y, z) has some known (exact) form by which we 
can check our numerical results. 

Rather than begin with this somewhat complicated equation, we begin by 
solving a slightly simpler-looking equation in two and three dimensions, in which 
the nonlinearity is of a quadratic form, as in 

V'^u{x, y, z) + ai?(x, y, z) = /(x, y, z), (21) 



Multigrid methods on domains containing holes 



7 



where tr (g M) is some tunable parameter which we typically set to ±1 (although there 
are good reasons to prefer —1 over +1; see footnote regarding smoothing operations, 
on the next page). We chose this equation not because of any particular physical 
relevance, but simply because it was the equation used in the version of Choptuik's 
2D multigrid code |35| . with which we had prior experience. 

Before proceeding to 3D calculations, we will describe our implementation for 
solving a 2D version of H21|l . Due to memory limitations (and array-indexing 
limitations in many Fortran compilers) , we can run our (non-parallel) 3D code only at 
much lower resolutions than we are able to achieve with the 2D code. Since some the 
interesting features in the convergence studies appear only only at very high resolutions 
in our 2D results, we first present our 2D scheme and the results yielded by it. 

3.1. Two dimensions 

Our starting point is a 2D FAS multigrid solver by Choptuik ^^1, which solves the 
equation 

-g^u{x, y) + ^w(a;, y) + <7u^{x, y) = f{x, y), (22) 

on a domain ^l with coordinate ranges [0, 0] to [1, 1], and subject to Dirichlet conditions 
at the outer boundary: u{x,y)\dno — const., and we choose this constant to be zero. 
The fimction f{x,y) is chosen such that the solution is 

u{x, y) = s\ii{Tilxx) sin(7r/yj/), (23) 

where Ix and ly are integers (which we set to unity for the cases presented in this 
paper). We also choose cr = 1. Choptuik's form assumed a domain without holes, but 
we extend this to configurations with holes, and on these inner boundaries we also 
apply Dirichlet conditions, the values for which are discussed below. 

The code uses a hierarchy of so-called "vertex centered" grids. Each grid at 
multigrid level I is a square lattice having 2' + 1 grid points along each edge. The 
grids have uniform spacing hi — 2~' in both x and y directions, and the grid 
points are denoted with indices i and j in the x and y directions, respectively, e.g., 
u{ihi,jhi) ~ u'^j. 

To Choptuik's original code, we added added features needed to handle the 
excision region. We also performed a few minor optimizations. At each grid point, we 
also store an integer code which denotes whether the point is a boundary point, an 
excised point, or an interior (or active) point. The set of all excised points is known 
as the excision mask or mask, and for the cases considered in this paper, we define the 
mask to be those points which are a distance r < Tmask from the center of the grid, 
where rmask is some value of our own choosing. The outermost points of the mask are 
where the inner boundary conditions are applied; we refer to these points as "inner 
boundary points." In addition to the single-hole case, we will consider the case of two 
holes, for which the excision regions will be points within the radius Tmask from the 
center of each hole. 

Boundary Conditions 

We use Dirichlet boundary conditions at both the inner boundary (edge of the 
hole) and the outer boundary, and apply them on all grid levels. At the outer boundary. 
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Figure 1. Example of how the inner boundary is defined, showing points on 
a coarse grid and a fine grid. Inner boundary points are those points which 
are immediately interior to a circle of radius rmax- The large filled circles show 
normal interior grid points (i.e., non-excised, non-boundary points) on the coarse 
grid, and the large open circles show boundary points on the coarse grid. The 
small filled and open circles show fine grid interior points and boundary points, 
respectively. The small dots show excised points on the fine and/or coarse grids, 
as appropriate. (Only one quadrant of a full domain is shown in this picture 
for purposes of clarity, but our simulations arc usually performed with the circle 
being centered on the grid.) 

we simply set = 0. At the inner boundary, we use the values of the exact solution, 
i.e., we set — u{x, y). 

The inner boundary is given by the outer edge of the points comprising the 
excision mask; an example is shown in Figure ^ Thus, in some sense these boundary 
points are not truly "excised" since they have data on them. This has the consequence 
that the "size" of the excision region (i.e., the area of the convex hull of the data 
points comprising the mask) on finer grids is always equal to or greater than the size 
of the excision region on coarser grids. This choice of the inner boundary, as being 
those points just inside some level surface instead of just outside it, is somewhat at 
variance from other excision schemes used in numerical relativity (e.g., |36p. Our 
choice is not fundamentally based on physical or mathematical principles, and is 
thus somewhat arbitrary. However we find that this definition of the inner boundary 
appears to be a key to our results of second order convergence: if we instead define 
inner boundary points to be those just beyond r = rmax, then we do not obtain second 
order convergence. 

Smoothing Operations 

We use a typical "red-black" Gauss-Seidel Newton iteration to smooth the error 
on each grid, in which we loop over all 'interior' points (i.e. non-excised, non-boundary 
points) and apply the following two equations % 

X We note an example of the care that must be taken in any solution of a discretized problem. 
Convergence of the solution may be impossible to achieve; this may signal significant analytical 
defects in the formulation. If ii > 0, then cr = — 1 guarantees that the Jacobian 2aUij — 4//i^ is 
nonzero (its inverse is nonsingular in Eq. 1251 '). Clearly, there are in principle situations with cru > 
where the Jacobian can vanish. Further, because the equations we wish to solve are nonlinear, zeros 
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Figure 2. 1-D schematic of scheme for inner boundary and restriction scheme. 
The circles on the rightmost X's indicate that this is where the Dirichlet conditions 
are applied, i.e., there are data on these points for all quantities except C^u^ 
(which appears in 1151 ^. which we do not compute on the boundary. One could 
use these points in weighted restriction even in the case shown in the right panel. 
However, we choose never to use these boundary data in weighted restriction, and 
instead do a simple "copy" operation. 



+ - (24) 
"^T-^f-^-J^- (25) 



Restriction and Prolongation Operators 



The restriction operator J^'* we use is the so-called "half-weighted" average on 
normal interior points, in which coarse grid values (indexed by / and J for clarity) are 
a weighted average of the fine grid values over a nearby region of the physical domain: 

= II" ^" = + I [«'Vi,, + ^-1, + "',+1 + "t-i] ' (26) 

where i = 21 — 1 and j = 2J — 1. Along the outer boundary, we perform a simple 
copy operation, Uj'j = u'^j. We use H26II for all restriction operations on interior grid 
points with one exception: We only use weighted restriction if none of the fine grid 
points used in the restriction operator is on or inside the boundary of the excised 
region; otherwise we use a "simple injection" or "copy" operation. This is illustrated 
in Figure [3 In conjunction with our definition of the inner boundary location, this 
use of a "copy" operation near the inner boundary is the central insight for preserving 
second order accuracy near the excision region. 

For the prolongation operator /j^^ , we use simple bilinear interpolation. 



3.1.1. 2D Results 



One of the first things we notice when graphing preliminary numerical solutions is 
that the solution error e = u — ii^ seems to be largest and "most in need of smoothing" 

of the Jacobian depend on ii and may occur at isolated spatial points. Further, because 2aUij 
is similar on all grids, but —Xjyfi is resolution dependent, this problem can appear differently on 
different multigrid levels. It may be that a more sophisticated solver could integrate through those 
singular points, but we do not address that question here. Notice that because Ui^j is near unity in 
the case of "sine" data given by Eq. I23i . this problem does not arise for tr = 1 for "sine" data. 
Similar comments apply, of course, in the 3D case of equation 1211 . 
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Figure 3. A plot of the solution error e = u—u for a system physical parameters 
a = 1 and Ix = ly = 1- We use a central circular hole of radius r-mask = 0.129, 
where this number is chosen such that the excision masks on different grid levels 
do not match up with one another, i.e., to make the solution harder to obtain, 
in order to demonstrate the robustness of our multigrid scheme. We see that 
the error is largest and sharpest next to the circular hole in the center. Thus, 
concentrating our smoothing operations near the hole is likely to increase the 
efficiency of the solver. These data are obtained from a run with parameters 
imax = 8, Lmin = 2, 1 V-cycle and 2 pre- and 2 post-CGC smoothing sweeps. 



(i.e., having high-frequency components) in the immediate vicinity of the excision 
region, as shown in Figure 01 

Thus we may perform only a few smooths on the entire grid and apply extra 
smoothing runs in an "extra smoothing region" (ESR) around the excision region, 
such as that shown in in the left panel of Figure^ In our implementation, we choose 
the width of the ESR to be a number of grid points which is one less than the current 
multigrid level number (e.g., on a level 4 grid, the ESR has a width of 3 grid points), 
and we smooth over the ESR twice as often as over the rest of the domain. We can 
compare the results of using the extra smoothing region to the results of smoothing 
over the entire domain. Such a comparison is shown in the right panel of Figure^ In 
our multigrid implementation, we always do twice as many smooths in the ESR as in 
the rest of the interior. 

Two Holes 

One of our principal applications for the multigrid solver will be to solve the 
constraint equations for binary black hole spacetimes. Thus we do a check to make 
sure that, at least for the nonlinear Poisson equation, our general method properly 
handles domains with multiple holes. One such case is shown in Figure 5. 

There is a contrast between the relative importance of pre-CGC smoothing versus 
post-CGC smoothing. Pre-CGC smoothing prepares the data in such a way that 
restriction will not produce aliasing errors, whereas post-CGC smoothing in large part 
corrects for errors brought in via the use of (linear) interpolation. It seems plausible 
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Figure 4. Left panel: Schematic of a level = igrid (2'"^^°' + 1 grid points per 
side) showing excised points (X's), normal interior points (N's), boundary points 
(B's) and the Extra Smoothing Region (ESR, E's). In this case, the excision 
region has radius rm^sk = 0.129. Right panel: Convergence behavior of Eqs. 
l22i - I^Sl with and without the use of the ESR. Here we show the L2 norm of 
the solution error e = u — u. We assign the width of ESR to be level — 1 grid 
points on either side of the hole, and we smooth twice as often over the ESR as 
over the normal interior points. Since the error is concentrated near the hole, we 
see that using 2 pre- and post-CGC smooths with the ESR (i.e., using 4 smooths 
in the ESR, and 2 smooths elsewhere) offers a substantial improvement over the 
use of 3 pre- and post-CGC smooths over the whole domain. We see that at high 
resolutions (level = 9 through 11) the error begins to rise again. To obtain results 
in keeping with second order convergence at resolutions up to level = 11, more 
work is required. The use of 3 pre- and post-CGC smooths with the ESR (i.e., 
6 smooths in the ESR, and 3 elsewhere) produces results comparable to those 
obtained by using 6 pre- and post-CGC smooths over the whole domain, however 
the former requires a fraction of the computational cost of the latter. 



that the former case (pre-CGC) would require at least as much smoothing as the 
latter (post-CGC); we see evidence for this in the two-hole solutions, such as those 
used to produce Figure|Sl Again we are able to achieve second order convergence, for 
instance with 1 V-cycle, 12 pre-CGC sweeps and no post-CGC sweeps, and with more 
effort: 5 V-cycles, 20 pre- and post-CGC sweeps, shown for comparison as a "perfect" 
multigrid result. 

3.2. Three dimensions 

For solving equations in 3D, we proceeded by steps, performing a few test cases 
until arriving at a code which solves an equation which is similar to the Hamiltonian 
constraint in 3-1-1 general relativity. 

The first test case is similar to the 2D case discussed above, i.e., we choose 
f{x, y, z) such that the exact solution to H21(l is 

u{x, y, z) — sin(7r^a;x) sin(7rZj,y) sin(7r^zz) (27) 

where , ly and 1^ are integers which we set to unity. We excised a sphere of radius 
''mask — 0.129, which is a value chosen simply to ensure that the size of the mask is 
different on different grid levels, i.e. to make the test a bit more difficult than if we 
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Figure 5. Comparison of pre-CGC smoothing vs. post-CGC smoothing, for 
a domain with two circular holes of radius rinask=0.1, separated by a distance 
d = 0.2, again for Eqs. 1221 - 1251 . Since a pre-CGC smooth requires the same 
number of operations as a post-CGC smooth, runs with 1 V-cycle, 3 pre- and 
3 post-CGC smooths are as costly as runs with 1 V-cycle, 6 pre- and post- 
CGC smooths. However the latter parameters yield a somewhat lower error at 
all resolutions. In another case, we compare 6 pre- and 6-post CGC smooths 
with 12 pre- and post-CGC smooths, with the similar result that 12 pre- and 
post-CGC smooths yield lower error. We note that the general shapes of these 
graphs seem to alternate between 'turning down' and 'turning up', with the cases 
involving more computational cost producing results closer to the ideal line. This 
suggests that, for a given range of resolutions, one is able to obtain results closer 
and closer to second order convergence as one increases the computational cost; an 
'extreme' case of 5 V-cycles, with 20 pre- and 20 post-CGC smooths is shown to 
illustrate this. It is important to note that post-CGC smoothing is not negligible 
in all cases: One needs post-CGC smooths particularly for runs involving multiple 
V-cycles in order to ensure convergence, e.g., a 3 V-cycle run with 12 pre- and 
post-CGC smooths immediately begins to 'blow up' as the resolution is increased. 
Adding post-CGC smooths to such a run considerably improves the convergence. 



had chosen some multiple of the mesh spacing. For the test case involving a sinusoidal 
solution (|27|l . the results are very similar to those of the 2D case. Having verified that 
the multigrid solver could adequately handle such a system and achieve second order 
accuracy, we proceeded to the next step. 

The next test case involves an equation H2()(l which corresponds to the solution of 
Hamiltonian constraint equation for conformally flat spacetimes svstem;37]: 

V^u - K^u^ + = (28) 
u' 

where and are arbitrary positive real constants related to the expansion rate 
of the 3-space, and r is the usual radial coordinate. We adjust the matter source term 
f{x, y, z) such that the continuum solution is 

u{x,y,z) = l + 2M/r. (29) 

Here r is the usual radial coordinate, and we choose M — 1. The excision region 
contains r = and the function values in this region are not calculated. On the 
inner boundary (the boundary of the excision region) , we use values of the continuum 
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solution as Dirichlet conditions on only the finest grid, and on coarser grids, we use 
data copied or extrapolated from finer grid data. One can simply copy values from 
the fine grid to the coarse grid wherever the boundary points on the fine and coarse 
grids coincide. Wherever they do not coincide, i.e., wherever the coarse grid boundary 
points are 'interior' to the fine grid boundary points, we find it quite adequate to use 
quadratic extrapolation (in one dimension) from the neighboring (non-excised) fine 
grid points to obtain the value at the location of a coarse grid boundary point. The 
results we obtain appear to be rather insensitive to the form of the extrapolation used 
(e.g., the shape of the stencil). The use of extrapolated fine-grid data is an important 
difference from the way inner boundary values are obtained in the 2D code, and does 
constitute a great improvement over the generic applicability of our method. 

For the outer boundary in this case, we can either apply Dirichlet conditions using 
the continuum solution u{x,y,z), or we can use a mixed boundary condition called 
the Robin condition. The Robin condition requires 

|-[r(u -1)]=0, (30) 

on the boundary. The Robin condition is the standard condition used in relativity 
when a solution with the form H29|) is desired (see, e.g., [HZl)- However M is not a 
known constant until the solution is known, so the condition which the solution (|29|l 
satisfies, namely (l3Uf) . is implemented as the boundary condition. 

There are several ways to implement the Robin condition. Rather than taking 
derivatives in the radial direction as required by l|5n|l . we follow Alcubierre |SH1 and 
instead take derivatives only in directions normal to the faces of our cubical domain. 
This approach is much simpler to implement than using radial derivatives. We 
implement this using first order difference stencils; for the -l-a; face of the domain, 
at location Xi, the stencil is simply 

Ui,j,k =^ 1 + {ui-i,j,k - l) ^'~^-^'''\ (31) 

where rij^k is the radial distance to lattice location {i,j, k). We have also implemented 
alternative stencils: second order stencils in the normal direction, and a first order 
stencil in the radial direction. 



Smoothing Operations 



We use a simple "red-black" Gauss-Seidel Newton iteration to smooth the error. 



tgs — h ^ {ui+ij^k + Wi-ij-,fc + Uij+i^k + Uij-i^k 

+Ui,j,k+1 + Uij^k-1 ~ GUij^k) 



K'ul,^+A'u-l,~f,,,,k (32) 



Ui 1 k = Ui i k : — 5 (33) 



Restriction and Prolongation Operators 



For the restriction operator if/^ , we use a weighted average which is defined via 
a 3 X 3 X 3 stencil, with a weight of 1/8 on the central point, 1/16 on the center of 
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each face, 1/32 on the center of each edge and 1/64 on each corner, i.e.: PO] 

~2h _ j2h ~h _ 1 r-/i 

+ 2 (""i+l + ^i-l.j,k + j- + l,/c 

~^''^i+l,j,k+l + ''^i'-l j,fe+l + ^i+l,j,k-l + ""i-lj.fc-l 
+^i+l,i+l,fc + ""i-lj + l.fe + '^i+l,i-l,A; + ^i-l,i~l,fc) 

+ g (""i+lj + l.fc+l + ''^i-lj + l.fc+l + ^i+l,j-l,fe+l 

+'^f-lj-l,fe+l + ^i+l,j + l,A;-l + ""f-lj + l.fe-l 

+ + (34) 

As in the 2D case, if any of the find grid points included in the weighting are inner 
boundary points, then we use a simple copy instead of weighted restriction 
For the prolongation operator /j;^ , we use simple trilinear interpolation. 

3.3. 3D Results 

A graph of the convergence behavior for the solution of the conformally flat background 
Hamiltonian constraint, Eq. H28(l . is given in Figure The left panel shows that we 
obtain lines which run parallel to the line for second order behavior at high resolutions, 
for both a Dirichlet outer boundary condition and the "first order perpendicular" 
implementation of the Robin condition suggested by Alcubierre |38| . In fact, the two 
graphs we obtain lie on top of each other. The right panel shows the convergent 
solution error of our solutions (i.e. the difference between our solutions and the exact, 
continuum solution) 

4. Toward More General Applications 

We have made use of a known exact solution in two key elements of this paper: (1) 
supplying values for Dirichlet conditions on the inner boundary, and (2) calculating 
the solution error e = u — vP' for measuring the accuracy and convergence of the code. 
Since we wish to use the multigrid solver for situations in which an exact solution is 
not known across the whole domain (such as in a solution for the conformal factor 
in a general binary black hole spacetime), we describe here the modifications to the 
previous discussion in the absence of an exact solution. 

In the 2D code we used the continuum solution to supply Dirichlet conditions on 
the inner boundaries of all multigrid levels, however for the 3D code we only used the 
continuum solution on the finest grid, and then extrapolated or copied data from finer 
grids to coarser grids. (There was nothing about the 2D case that made extrapolation 
impossible; rather it was simply a later feature which was added to the more advanced 
3D code.) Although we were able to employ this extrapolation effectively in the cases 
we tried — using a very simple, almost arbitrary, ID extrapolation method — it may 
be that such a method would be unstable for certain classes of equations. Although 
we see no evidence for this, and given the generality of the multigrid method one may 
expect much of what is found in this paper to apply to other systems of interest. We 
cannot offer complete assurance that extrapolation will work for all elliptic systems. 



J 




Level 



Figure 6. Convergence results for 3D solutions to 1281 of the form u = l + 2M/r, 
for runs in which = = 1, l^^in = 3, and rnjask = 1-29, for a domain 
[—5, —5, —5] to [5,5,5]. Left panel: A logarithmic plot of the L2 norm of 
the solution error e = u — , showing a comparison between outer boundary 
conditions. Using the "first order perpendicular" (FOP) implementation 1311 of 
the Robin boundary condition I30i 'l. we obtain convergence results which lie on 
top of those obtained using a Dirichlet outer boundary condition (where the values 
of the continuum solution are supplied at the outer boundary). These results 
also run parallel to the line for perfect second order behavior. Right Panel: 
A logarithmic plot of e itself, at the end of each V-cycle in the Full Multigrid 
Algorithm. These results correspond to the "IV, 3 Pre, 3 Post, 1st order Perp" 
case shown in the left panel. Here we show a ID slice along the x-axis, and we have 
divided coarser grid values by appropriate powers of four in order to make the 
comparison. We see second order convergence over the interior of the domain and, 
importantly, in the vicinity of the excision region. Near the outer boundary, the 
magnitude of the error is roughly second order (thus it does not noticeably affect 
the graph shown in the left panel), however its shape is resolution-dependent. 
This feature may arise from the use of the FOP condition. 

The assumption that an exact value is known at points on the finest grid may 
not be appropriate for many of the solutions of interest to numerical relativists. One 
may encounter systems with Neumann or Robin conditions at the inner boundary. 
Although we have not considered such cases explicitly, we suggest straightforward 
differencing similar to the treatment of the Robin conditions we reported above. 

When one cannot measure the accuracy of the code by calculating the solution 
error, one can still gain a measure of the convergence of the code by monitoring the 
relative truncation error t^^ . One should disregard values at points adjacent to inner 
boundary points, as these tend to be poorly defined and exhibit blow-up behavior, 
but all other points are eligible for comparison. A plot of t^^ at various grid levels for 
a 3D solution is shown in Figure [7| demonstrating the strict second order behaviour 
of the solution on two different domains. 

5. Conclusions 

We can offer five tips for those who wish to implement 2D or 3D multigrid methods 
for domains with holes. 

• 1) The first and most important tip is really twofold: 
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Figure 7. A ID slice through the center of a 3D domain, for the truncation error 
of a solution to Eq. I28i with Robin outer boundary conditions, seeking a solution 
of the form u = l+2Af/r. The run parameters , and M are set to unity, and 
we perform 1 V-cycle, 3 pre-CGC smooths and 3 post-CGC smooths. Left panel: 
On the domain [—5, —5, —5] to [5,5,5], with r^^^j. = 1.29. Right panel: On 
the domain [-0.5,-0.5,-0.5] to [0.5,0.5,0.5], with r^ask = 0.129. (This smaller 
domain is intended to provide a stronger test of the algorithm than the larger 
domain.) The notation in the figure legends denotes which grid levels were 
used to calculated t^^ for a given set of points. We see that, on both domains, the 
0{h'^) error function is independent of resolution, in accordance with the idea of 
Richardson extrapolation. Since the continuum solution is u{x, y,z) = 1 + 2M/r, 
we expect the error function to tend toward large values as r — > 0. In these graphs, 
we have cut out the values of t'^^ immediately adjacent to the inner boundary. 



— apply inner boundary conditions on points immediately interior to (rather 
than exterior to, as is often done) some spatial surface such as a circle or 
sphere of a given radius. This has the effect that the extent of excision 
region is smaller on coarser grids than on finer grids. 

— when performing restriction operations, use weighted restrictions on all 
interior points, except in cases where the weighted operator would include 
points where the inner boundary conditions are applied. For these latter 
cases, use a simple copy operation. 

• 2) Concentrate smoothing operations where they are most needed. We have done 
this by arbitrarily defining an "extra smoothing region" around any hole, and 
performing twice as many smooths in this region as in the rest of the domain. 
One can imagine more sophisticated schemes, which look at derivatives of the 
solution or the local truncation error as indicators of where the extra smoothing 
should be performed and how much extra work should be performed there. 

• 3) It may be the case that pre-CGC smoothing is much more effective than post- 
CGC smoothing for a given problem, particularly if only one V-cycle is performed. 
Concentrating the bulk of the smoothing operations into the pre-CGC smooths 
results in a faster route to the solution at a desired accuracy. 

• 4) This repeats the advice of Alcubierre (2S| (Radial) Robin boundary conditions 
can be adequately mimicked by applying them only in the normal directions, thus 
simplifying the computer code and reducing the computation effort slightly. The 
Robin conditions should be applied at all grid levels; i.e. it is not the case that 
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one should apply Robin conditions on the finest grid only and then use those 
values as Dirichlet conditions on coarser grids. 

• 5) This tip is intended for use in situations where an exact solution is not available; 
we have less experience with these situations. However, we find that a naive 
quadratic extrapolation from nearby points is adequate to provide values for 
Dirichlet boundary conditions on the edge of the hole. Further (more technically), 
when monitoring the relative truncation error r^'', points immediately adjacent 
to inner boundary points have erratic behavior, and should not be considered in 
judging the behavior of the solution. 

We clearly have some further work ahead in order to consider more general 
systems in which our inner boundary data is not so easily obtained, and we also need to 
consider systems more complex than the Hamiltonian constraint for a conformally flat 
background, as treated in this paper. It is thus our hope that the techniques described 
here will also be of relevance for more interesting systems of elliptic equations. We 
are developing a solver for the full set of constraint equations (testing it using the 
Kerr-Schild-type binary black hole data described in '39' ) , and see promising results 
using precisely the techniques given in this paper. Further improvements in the speed 
of multigrid scheme, such as the use parallelism and adaptivity |4U[ 1411 will also 
likely be necessary for its practical application in a constrained evolution code. 

The 2D and 3D codes described here are available for use by the community, 
by download from littp://wwwrel. ph. utexas.edu/~shawley/mg_ex. html. We 
encourage other users of our approach. However, the codes are "as is" , and we cannot 
maintain them, nor correct problems that may arise from their use. 
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